R Notebook

This is an R Markdown Notebook to study boundary layer meteorology.

Create the dataset - - From Great Plains Field Prpgram at O’Neil Nebraska

TatDpoint025 <- c(25.54,24.84,25.77,29.42,33.25,35.25,34.84,32.63,30.07,28.42,27.09,26.09,25.30)
TatDpoint05 <- c(26.00,25.30,25.35,27.36,30.32,32.63,33.20,32.05,30.20,28.74,27.50,26.60,25.83)
TatDpoint10 <- c(26.42,25.84,25.42,25.98,27.62,29.52,30.62,30.62,29.91,28.84,27.84,27.06,26.40)
TatDpoint20 <- c(26.32,25.97,25.56,25.39,25.57,26.11,26.88,27.41,27.68,27.57,27.22,26.87,26.53)
TatDpoint40 <- c(24.50,24.48,24.36,24.34,24.27,24.24,24.26,24.32,24.47,24.64,24.73,24.78,24.84)
datetime <- c('2015-08-31 04:35','2015-08-31 06:35','2015-08-31 08:35','2015-08-31 10:35','2015-08-31 12:35','2015-08-31 14:35','2015-08-31 16:35','2015-08-31 18:35','2015-08-31 20:35','2015-08-31 22:35','2015-09-01 00:35','2015-09-01 02:35','2015-09-01 04:35')   
blaygpfp <- data.frame(TatDpoint025,TatDpoint05,TatDpoint10,TatDpoint20,TatDpoint40,datetime)

check the data

#convert the date format
blaygpfp$datetime<- as.POSIXct(blaygpfp$datetime, "%Y-%m-%d %H:%M",tz="US/Pacific")

unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'Warning in as.POSIXlt.POSIXct(x, tz) :
  unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

#TODO Take the annnoying date fotmat error out
str(blaygpfp)

'data.frame':   13 obs. of  6 variables:
 $ TatDpoint025: num  25.5 24.8 25.8 29.4 33.2 ...
 $ TatDpoint05 : num  26 25.3 25.4 27.4 30.3 ...
 $ TatDpoint10 : num  26.4 25.8 25.4 26 27.6 ...
 $ TatDpoint20 : num  26.3 26 25.6 25.4 25.6 ...
 $ TatDpoint40 : num  24.5 24.5 24.4 24.3 24.3 ...
 $ datetime    : POSIXct, format: 

unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

"2015-08-31 04:35:00" "2015-08-31 06:35:00" "2015-08-31 08:35:00" ...

Plot the data

plot(blaygpfp)

library(ggplot2)
library(dplyr)
library(reshape)

Attaching package: ‘reshape’

The following object is masked from ‘package:dplyr’:

    rename

Plot Temperature waves as function of time and depth.

blaygpfp %>% ggplot() +
geom_point(aes(x=datetime,y=TatDpoint025,colour = as.factor(0.025), shape=as.factor(0))) +
  stat_smooth(aes(x = datetime, y = TatDpoint025), se = FALSE) +
geom_point(aes(x=datetime,y=TatDpoint05,colour = as.factor(0.05))) +
  stat_smooth(aes(x = datetime, y = TatDpoint05), se = FALSE) +
geom_point(aes(x=datetime,y=TatDpoint10,colour = as.factor(0.10))) +
  stat_smooth(aes(x = datetime, y = TatDpoint10), se = FALSE) +
geom_point(aes(x=datetime,y=TatDpoint20,colour = as.factor(0.20))) +
stat_smooth(aes(x = datetime, y = TatDpoint20), se = FALSE) +
geom_point(aes(x=datetime,y=TatDpoint40,colour = as.factor(0.40))) +
stat_smooth(aes(x = datetime, y = TatDpoint40), se = FALSE) +
labs(colour = "Depth")  +
xlab("Date") + ylab("Temperature")

# check if line plot is better

With markers

#normalize temperature
blaygpfp %>% ggplot() +
geom_point(aes(x=datetime,y=TatDpoint025, shape=as.factor(0.025))) +
  geom_line(aes(x = datetime, y = TatDpoint025)) +
geom_point(aes(x=datetime,y=TatDpoint05,shape = as.factor(0.05))) +
  geom_line(aes(x = datetime, y = TatDpoint05)) +
geom_point(aes(x=datetime,y=TatDpoint10,shape = as.factor(0.10))) +
  geom_line(aes(x = datetime, y = TatDpoint10),) +
geom_point(aes(x=datetime,y=TatDpoint20,shape = as.factor(0.20))) +
geom_line(aes(x = datetime, y = TatDpoint20)) +
geom_point(aes(x=datetime,y=TatDpoint40,shape = as.factor(0.40))) +
geom_line(aes(x = datetime, y = TatDpoint40)) +
labs(colour = "Depth")  +
xlab("Date") + ylab("Temperature")+
ggtitle("Temperatures waves - function of time and depth")  

# check if line plot is better

Housekeeping

# reshape it
blaygpfp_melt <- melt(blaygpfp,id="datetime")

Warning in as.POSIXlt.POSIXct(x, tz) :
  unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

# set the depth
blaygpfp_melt$depth <- 0

Warning in as.POSIXlt.POSIXct(x, tz) :
  unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

blaygpfp_melt[blaygpfp_melt$variable=='TatDpoint05',]$depth = 0.05

Warning in as.POSIXlt.POSIXct(x, tz) :
  unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

blaygpfp_melt[blaygpfp_melt$variable=='TatDpoint10',]$depth = 0.10

Warning in as.POSIXlt.POSIXct(x, tz) :
  unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

blaygpfp_melt[blaygpfp_melt$variable=='TatDpoint20',]$depth = 0.20

Warning in as.POSIXlt.POSIXct(x, tz) :
  unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

blaygpfp_melt[blaygpfp_melt$variable=='TatDpoint40',]$depth = 0.40

Warning in as.POSIXlt.POSIXct(x, tz) :
  unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

blaygpfp_melt[blaygpfp_melt$variable=='TatDpoint025',]$depth = 0.025

Warning in as.POSIXlt.POSIXct(x, tz) :
  unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

library(lubridate)

Attaching package: ‘lubridate’

The following object is masked from ‘package:reshape’:

    stamp

The following object is masked from ‘package:base’:

    date

#In check_tzones(e1, e2) : 'tzone' attributes are inconsistent
blaygpfp_melt$datetime <- ymd_hms(blaygpfp_melt$datetime)

unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

#Do increment of the time intervals
blaygpfp_melt %>% filter(datetime == ymd_hms("2015-09-01 04:35:00")) 

# set the time intervals
 blaygpfp_melt$interval <- 0
blaygpfp_melt[datetime == ymd_hms("2015-08-31 06:35:00"), ]$interval<- 2
blaygpfp_melt[datetime == ymd_hms("2015-08-31 08:35:00"), ]$interval<- 4
blaygpfp_melt[datetime == ymd_hms("2015-08-31 10:35:00"), ]$interval<- 6
blaygpfp_melt[datetime == ymd_hms("2015-08-31 12:35:00"), ]$interval<- 8
blaygpfp_melt[datetime == ymd_hms("2015-08-31 14:35:00"), ]$interval<- 10
blaygpfp_melt[datetime == ymd_hms("2015-08-31 16:35:00"), ]$interval<- 12
blaygpfp_melt[datetime == ymd_hms("2015-08-31 18:35:00"), ]$interval<- 14
blaygpfp_melt[datetime == ymd_hms("2015-08-31 20:35:00"), ]$interval<- 16
blaygpfp_melt[datetime == ymd_hms("2015-08-31 22:35:00"), ]$interval<- 18
blaygpfp_melt[datetime == ymd_hms("2015-09-01 00:35:00"), ]$interval<- 20
blaygpfp_melt[datetime == ymd_hms("2015-09-01 02:35:00"), ]$interval<- 22
blaygpfp_melt[datetime == ymd_hms("2015-09-01 04:35:00"), ]$interval<- 24
#'data.frame':  65 obs. of  3 variables:
# $ datetime: POSIXct, format: "2015-08-31 04:35:00" "2015-08-31 06:35:00"  #$ variable: Factor w/ 5 levels "TatDpoint025",..: 1 1 1 1 1 1 1 1 1 1 ...
# $ value   : num  25.5 24.8 25.8 29.4 33.2 ...

Plot Time waves as function of depth and temperature

#normalize temperature
blaygpfp_melt$scaleT <- scale(blaygpfp_melt$value)  %>% as.vector()
# We want for hours 0, 4,8,12,16,20
blaygpfp_melt  %>% filter(interval %in% c(0, 4,8,12,16,20)) %>%  ggplot() +
geom_point(aes(x=value,y=depth,color=as.factor(interval),shape=as.factor(interval)))+
   
#geom_line(aes(x=value,y=depth,color=as.factor(interval),shape=as.factor(interval)))+
stat_smooth(aes(x=value,y=depth,color=as.factor(interval),shape=as.factor(interval)))+  
xlab("Temperature(deg C)") + ylab("Depth")+
ggtitle("Temperatures at various times - function of z and T")  

Ignoring unknown aesthetics: shape

# check if line plot is better

plot for depth

library(ggplot2)
amp <- c(5,3.6,2.5,1.5,.5)
depth <- c(0.025,0.05,.1,.2,.4)
plotampdepthdf <- data.frame(depth,amp)
lm_eqn <- function(df){
    m <- lm(amp ~ depth, df);
    eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
         list(a = format(coef(m)[1], digits = 2), 
              b = format(coef(m)[2], digits = 2), 
             r2 = format(summary(m)$r.squared, digits = 3)))
    as.character(as.expression(eq));                 
}
p <- ggplot(data = plotampdepthdf, aes(x = depth, y = amp)) +
            geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
            geom_point()
p1 <- p + geom_text(x = 0.2, y = 3, label = lm_eqn(plotampdepthdf), parse = TRUE)
p1

With log Amp

library(ggplot2)
amp <- c(5,3.6,2.5,1.5,.5)
depth <- c(0.025,0.05,.1,.2,.4)
amp <- log(amp)
plotampdepthdf <- data.frame(depth,amp)
lm_eqn <- function(df){
    m <- lm(amp ~ depth, df);
    eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
         list(a = format(coef(m)[1], digits = 2), 
              b = format(coef(m)[2], digits = 2), 
             r2 = format(summary(m)$r.squared, digits = 3)))
    as.character(as.expression(eq));                 
}
p <- ggplot(data = plotampdepthdf, aes(x = depth, y = amp)) +
            geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
            geom_point()
p1 <- p + geom_text(x = 0.2, y = .5, label = lm_eqn(plotampdepthdf), parse = TRUE)+
    xlab("Depth") + ylab("log(amplitude)") + ggtitle("Log of Amplitude")
p1

Lets try a polynomial fit

lm(y ~ poly(x, 3, raw=TRUE))

amp <- c(5,3.6,2.5,1.5,.5)
depth <- c(0.025,0.05,.1,.2,.4)
plotampdepthdf <- data.frame(depth,amp)
poly.model <-lm(amp ~ poly(depth,3,raw=TRUE), df)
summary(poly.model)

Call:
lm(formula = amp ~ poly(depth, 3, raw = TRUE), data = df)

Residuals:
        1         2         3         4         5 
 0.111083 -0.208280  0.121497 -0.026035  0.001736 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)  
(Intercept)                    6.2571     0.4958  12.620   0.0503 .
poly(depth, 3, raw = TRUE)1  -60.9184    13.9633  -4.363   0.1434  
poly(depth, 3, raw = TRUE)2  256.3255    89.8197   2.854   0.2146  
poly(depth, 3, raw = TRUE)3 -350.0562   146.3810  -2.391   0.2521  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2668 on 1 degrees of freedom
Multiple R-squared:  0.9943,    Adjusted R-squared:  0.977 
F-statistic: 57.69 on 3 and 1 DF,  p-value: 0.09641

lm_eqn <- function(df){
    m <- lm(amp ~ poly(depth,3,raw=TRUE), df);
    eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
         list(a = format(coef(m)[1], digits = 2), 
              b = format(coef(m)[2], digits = 2), 
             r2 = format(summary(m)$r.squared, digits = 3)))
    as.character(as.expression(eq));                 
}
p <- ggplot(data = plotampdepthdf, aes(x = depth, y = amp)) +
            geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
            geom_point()
p1 <- p + geom_text(x = 0.2, y = 3, label = lm_eqn(plotampdepthdf), parse = TRUE) +
    xlab("Depth") + ylab("Amplitude") + ggtitle("Polynomial fit(Order 3)")
p1

Lets try to plot the predicted

plm_eqn <- function(df){
    m <- lm(amp ~ poly(depth,3,raw=TRUE), df);
    eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
         list(a = format(coef(m)[1], digits = 2), 
              b = format(coef(m)[2], digits = 2), 
             r2 = format(summary(m)$r.squared, digits = 3)))
    as.character(as.expression(eq));                 
}
predictedamps <- predict(poly.model,list(depth=plotampdepthdf$depth))
p <- ggplot(data = plotampdepthdf, aes(x = depth, y = amp)) +
            geom_point()
p1 <- p + geom_line(aes(x=depth,y=predictedamps)) +
     geom_text(x = 0.2, y = 3, label = plm_eqn(plotampdepthdf), parse = TRUE) +
     xlab("Depth") + ylab("Amplitude") + ggtitle("Predicted with actual")
p1

TODO - plot the tl vs depth

library(ggplot2)
lag <- c(5,3.6,2.5,1.5,.5)
depth <- c(0.025,0.05,.1,.2,.4)
plotampdepthdf <- data.frame(depth,amp)
lm_eqn <- function(df){
    m <- lm(amp ~ depth, df);
    eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
         list(a = format(coef(m)[1], digits = 2), 
              b = format(coef(m)[2], digits = 2), 
             r2 = format(summary(m)$r.squared, digits = 3)))
    as.character(as.expression(eq));                 
}
p <- ggplot(data = plotampdepthdf, aes(x = depth, y = amp)) +
            geom_smooth(method = "lm", se=FALSE, color="black", formula = y ~ x) +
            geom_point()
p1 <- p + geom_text(x = 0.2, y = 3, label = lm_eqn(plotampdepthdf), parse = TRUE)
p1

We need the time at which the maxima occurs, so do the lazy way to get the maxima by installing plotly to get visual cues.

blaygpfp %>% ggplot() +
  stat_smooth(aes(x = datetime, y = TatDpoint025,colour = as.factor(0.025)), se = FALSE) +
  stat_smooth(aes(x = datetime, y = TatDpoint05,colour = as.factor(0.05)), se = FALSE) +
  stat_smooth(aes(x = datetime, y = TatDpoint10,colour = as.factor(0.10)), se = FALSE) +
stat_smooth(aes(x = datetime, y = TatDpoint20,colour = as.factor(0.20)), se = FALSE) +
stat_smooth(aes(x = datetime, y = TatDpoint40,colour = as.factor(0.40)), se = FALSE) +
labs(colour = "Depth(m")  +
xlab("Date") + ylab("Temperature (deg C)")+
ggtitle("Soil Temperatures at various depths GPFP (O'Neil Nebraska) ")  

# check if line plot is better
ggplotly()

`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

Redo as line , but no smoothing

blaygpfp %>% ggplot() +
geom_line(aes(x=datetime,y=TatDpoint025,colour = as.factor(0.025))) +
geom_line(aes(x=datetime,y=TatDpoint05,colour = as.factor(0.05))) +
geom_line(aes(x=datetime,y=TatDpoint10,colour = as.factor(0.10))) +
geom_line(aes(x=datetime,y=TatDpoint20,colour = as.factor(0.20))) +
  geom_line(aes(x=datetime,y=TatDpoint40,colour = as.factor(0.40))) +
labs(colour = "Depth(m")  +
xlab("Date") + ylab("Temperature (deg C)")+
ggtitle("Soil Temperatures at various depths GPFP (O'Neil Nebraska) ")  

# check if line plot is better
ggplotly()

unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'unknown timezone 'zone/tz/2017c.1.0/zoneinfo/America/Los_Angeles'

Depth and time of the temperature

maxtempat <- c(14,16,16,20,4)
depth <- c(0.025,0.05,.1,.2,.4)
plothourdepthdf <- data.frame(depth,maxtempat)
p <- ggplot(data = plothourdepthdf, aes(x = depth, y = maxtempat)) +
            geom_point()
p

Very odd plot

try with hour start at 0

TODO - use - blaygpfp_melt df

Vertical soil temperatures at 0435, 0835, 1235,1635,2035 and 0035

library(lubridate)

fix the timestamp